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We show by atomistic simulations that, in the thermodynamic limit, the in-plane elastic moduli of 
graphene at finite temperature vanish with system size L as a power law L“’'“ with rju — 0.325, in 
agreement with the membrane theory. Our simulations clearly reveal the size and strain dependence 
of graphene’s elastic moduli, allowing comparison to experimental data. Although the recently 
measured difference of a factor 2 between the asymptotic value of the Young modulus for tensilely 
strained systems and the value from ab initio calculations remains unsolved, our results do explain 
the experimentally observed increase of more than a factor 2 for a tensile strain of only a few 
permille. We also discuss the scaling of the Poisson ratio, for which our simulations disagree with 
the predictions of the self-consistent screening approximation. 


Mechanical and structural properties of graphene form 
an intriguing and highly non-trivial aspect of its physics. 
The structure of a two-dimensional (2D) material em¬ 
bedded in a 3D space gives room to special features, 
related to large out-of-plane deformations, in particular 
thermal ripples m and static ripples and wrinkles m- 
A crucial difference with 3D (or strictly 2D) crystals is 
that, for graphene, out-of-plane atomic displacements h 
and in-plane displacements u have different wavevector 
dependence of the energy cost in the long wavelength 
limit g —)■ 0, namely oc and oc q respectively, the 
latter being the normal behavior for acoustic phonons. 
Hence, at finite temperature, the long wavelength out- 
of-plane fluctuations are much larger than the in-plane 
ones, so that at some small wavevector q the first anhar- 
monic coupling term of the form uh^ will dominate over 
the “normal” harmonic terms u^, with important con¬ 
sequences for the elastic behavior [SI 13 fTTHTi] . In par¬ 
ticular, the wavevector dependence of the anharmonic 
coupling strength leads to expect a power law behavior 
of the size dependence of the elastic properties. 

Contrary to the temperature dependence [13 (13 , so 
far the size dependence of the in-plane elastic moduli 
of graphene has been hardly studied nor measured m 
until recent experiments seem to indicate that such a 
size dependence does exist for graphene [I3II3- From 
indentation experiments on graphene drums with sizes 
of the order of 1 fj,m, the Young modulus Y was found 
to vary between 250 N/m and 700 N/m with increasing 
strain m, the latter value being much higher than the 
currently accepted value for flat sheets of ~ 340 N/m 
obtained in previous measurements HU or ab initio cal¬ 
culations [23 • 

Here we study the size dependence of the in-plane elas¬ 
tic moduli of graphene at room temperature r=300 K by 
means of atomistic Monte Carlo (MC) simulations based 
on the realistic interatomic potential LCBOPH [ 21 ], as 
used in previous works We obtain explicit ex¬ 

pressions for the size and strain dependence of graphene’s 
in-plane elastic moduli, providing a benchmark and tools 


for the analysis of experiments for systems of any size. 

Theoretically, the mentioned size dependence has been 
studied within the continuum elastic theory of thin plates 
and membranes, described by the Hamiltonian [22] : 

+ Xul^ + 2fiulp) ( 1 ) 


where r is the 2D position vector, k is the bending rigid¬ 
ity, A and /r are Lame coefficients, with /i the shear mod¬ 
ulus, and 

lla/3 — 2 T T (2) 


is the strain tensor. 

The harmonic approximation neglects the non-linear 
term, decoupling the bending and stretching modes. 
Then the correlation functions for out-of-plane displace¬ 
ments, Ho{q) = (|hqp)o, and for in-plane displacements, 
0 ( 9 ) = ('*^aq'*^/3q)oi can be derived by Gaussian inte¬ 
gration 


Ho{q) = 


keT 

nq'^ 
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and: 
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with Pa,a(q) = 9 q 9 / 3 / 9 ^- The average height fluctuation 
behaves as (h^)o = X]q(l^qP)o ~ implying instabil¬ 
ity of a membrane as a flat phase. 

Due to the large out-of-plane fluctuations, however, the 
harmonic behavior is not valid for small q and one has 
to keep the term in eq. Since H remains quadratic 
in u, these degrees of freedom can still be integrated out. 
This leads to a Hamiltonian in Fourier space which is a 
function of hq only [3: 
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where Y is the 2D Young modulus and i?(q, k, k') = (q x 
k)(qx k')/g^. The anharmonic, quartic term reduces the 
height fluctuations, stabilizing the flat phase, effectively 
described by a renormalized bending rigidity Kfi{q) ^ 
with positive rj. Hence, the height correlation H{q) 
for g —^ 0 has the same form as Ho{q) in eq. but 
with K replaced by Kji{q) [TT]. Likewise, D'^^[q) can be 
described by renormalized \fi{q), iiii{q) ^ q^'‘ in eq. 
with riu > 0. From rotational invariance it follows that 
rj and ??„ should satisfy the scaling relation = 2 — 2ri 

[23]. 

Within the self-consistent screening approximation 
(SCSA) [T3], the exponent was estimated as ry ~ 0.821 
m-: next-order corrections reduce it slightly to 77 ~ 0.789 
M- A renormalization group approach gives 77 = 0.849 
[24| and MC simulations for self-avoiding membranes 
77 ~ 0.72 [25]. With 77 > 0, (/i^) ^ is much smaller 

than (/i^)o ~ for large L, stabilizing the flat phase. 

Although it is a priori not obvious whether the mem¬ 
brane theory applies to an atomic-layer-thick 2D crys¬ 
tal like graphene, atomistic MC simulations confirm the 
scaling behavior of H{q) with 77 ~ 0.85|4j. The scaling of 
in-plane elastic moduli, however, has not yet been stud¬ 
ied nor confirmed for graphene. Contrary to kr which 
increases with increasing system size, making the mem¬ 
brane more resistant against bending, \r and pR de¬ 
crease with system size. Hence, if graphene follows the 
membrane theory, the in-plane elastic moduli vanish for 
large system sizes, an unthinkable situation for 3D crys¬ 
tals! 

For a 2D system, the 2D bulk modulus B, the uniaxial 
modulus Cii and Y are related to A and p as: 

B = X + p , Cii = B + p and Y = - — ( 6 ) 

B + p 

implying that H, Cn and Y scale as A and p. Another 
relevant quantity is the 2D Poisson ratio v: 


The SCSA predicts a universal, negative Poisson ra¬ 
tio iz = —1/3 for L —)■ 00 [T3], as later confirmed by 
MC simulation of self-avoiding membranes [2^. For 
graphene, however, so far only positive values were re¬ 
ported [v = 0.15 - 0.46) [ISII^Eg. 

In Fig [2 we show the in-plane correlation function 
D'^°‘{q) {a = x,y) calculated by NPT MC simulations 
at pressure P = 0 and T = 300 K with isotropic volume 
fluctuations for roughly square samples of A' = 37888 
atoms using periodic boundary conditions. Besides dis¬ 
placement moves we apply also collective wave moves for 
small q as in Ref. El For the calculation of D'^°‘{q) = 
(kgP) = ( 1 /A^)(IE*^ M,„exp(iqr,_o)| 2 ) with {r,,o} the 
ground state positions, the in-plane displacement field 
was scaled as Uia = svia — Tiap where s = ^A^jA scales 


the area A at T=300 K to the ground state area Aq of a 
flat sample. The behavior of D^^{q) for small q is consis- 
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FIG. 1. Correlation functions DZ°‘ {a = x,y) for in-plane 
displacements Uix (x) and Uiy (-I-). The scaling exponent is 
consistent with with rju = 2 — 2ri = 0.3, using 

77 ~ 0.85 Eg (dashed line). 

tent with a power law with exponent 77 ^ ~ 0.3, indicating 
that graphene follows the membrane theory also for in¬ 
plane correlations. 
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FIG. 2. Pressure as a function of size in NPT simulations 
(symbols) with (a) isotropic and (b) uniaxial size fluctuations 
for approximately square systems with N =24, 112, 336, 1008, 
4032, 12096 and 37888 atoms. The lines are best fits to eq. |11| 
The inset gives the equilibrium sizes Seq = s{P = 0) (symbols) 
as a function of Lg and the fit (solid line) according to the 
expression in Table jl] 

The size and strain dependence of the elastic properties 
can be computed simultaneously by NPT MC simulations 
for different sizes with isotropic area fluctuations at sev¬ 
eral pressures P. The resulting average area A gives the 
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equation of state (EOS), A{P) and thus also P{A) from 
which B can be calculated as: 


— 


2 ds 


( 8 ) 


where s = L/Lq is the relative linear system size with 
Lq = yjN /pq the ground state system size, po — 
0.3819 being the 2D ground state atomic density of 
graphene. To obtain both B and Cn, we also performed 
NPT simulations for uniaxial pressure P^, applying fluc¬ 
tuations of Lx in the ^-direction, while keeping Ly fixed. 
Then Cn follows from: 


C'li — —Lx 


dPx 

dPx 

dPx 

dLx 

ds 

p 

P ^^X 


(9) 


Sy — Seq 


where Sa = La/Lafi, with La^ {a = x,y) the ground 
state dimensions, and where Seq = s{P = 0) is the equi¬ 
librium size obtained from the isotropic NPT simulations 
at P = 0. The subscript “P’ in eq. indicates that 
Ly should be taken equal to Sy = s(P) resulting from 
isotropic NPT simulations at pressure P and that Px 
should be varied around P. However, since we verified 
that the dependence of dPx/dsx on Sy is very weak we 
adopted the last approximation in eq. which is exact 
for P = 0. 

The results are shown in Fig. The inset shows that 
the previously found negative thermal expansion m is 
also size dependent, but tending to a constant for large 
Lq. On the basis of Fig. with the slope dP/ds = 
2B/s tending to a constant for large s, we propose the 
phenomenological relation for B{s) 


Bis) = 


5 {Bgqf Seq p CDi^S Seq)) 

1 -I- L>(S - Seq) 


( 10 ) 


where Beq is the equilibrium value at P = 0. Substitution 
of eq. W into eq. [^and integration yields the EOS: 

-P(S) =-^ (—-C)ln(l+D(s- Seq)) - 2C(S - Seq) 


D 


( 11 ) 


Similarly, we write 


C'ii(sx) — 


(On 

,eq / ^eq + CD(Sx — Seq)^ 
1 p L)(^Sx ^eq) 


( 12 ) 


which substituted in eq. gives an equation for Px{sx) 
similar to eq. [TT] but with s, Beq, C and D replaced by 
Seq, Cii,eql‘^, 0/2 and D. This form allows the excellent 
fits shown in Fig. providing Beq and On,eg as a func¬ 
tion of Lq. In the left panels of Fig. (left panels) we 
show that both B and On vanish for large Lq, decreasing 
as a power law ^ P”’*", with rju — 0.325 (insets). The 
right panels give the corresponding results for Y and iz at 
P = 0, calculated using eqs. |^and[^ Note that, accord¬ 
ing to LCBOPII, the in-plane elastic moduli of graphene 



FIG. 3. Equilibrium bulk modulus B^q, uniaxial modulus 
On,eg, Young modulus Y^q and Poisson ratio iz as a function 
of system size Lq. The insets in log-log scale demonstrate the 
power law behavior. The solid lines are fits according to the 
expressions in Table |I] Dashed lines in the right panels are 
guides to the eye. 


at T = 0 K are P = 12.69 eV/A^ and p = 9.26 eV/AA 
yielding Y = 21.41 eV/A^ = 343 N/m and v = 0.156, in 
agreement with ab initio data |20] and with the small size 
limit in Fig. where Y ~ 314 N/m. By simulations at 
1 K for = 24 we verified that the remaining difference 
is due to temperature. 

Interestingly, the power law decrease of B, Cu and 
Ygq as a function of Lq sets in from Lq ~ 20 A, a 
value twice smaller than the Ginzburg critical value 
L* = 2tt16tt/{SY ksT) ~ 40 A (using k ~ l.leH [5]) 
expected from membrane theory m- The Poisson ratio 
1 / for small sizes is close to its bare value and increases 
up to 0.275 for larger Lq, against the SCSA prediction 
1 / = —1/3. Since the scaling of B and p is consistent with 
the SCSA, it is very unlikely that v will reach the value 
-1/3 for Lq — )• 00 , as the outcome of eq. [^only depends 
on the prefactors. 

We can also calculate Y as a function of tensile strain, 
using eqs. [T0| and [T^ with the best fit parameters. We 
should use the B{s) and Cii{sx) at equal pressure by 
solving Px{sx) = P{s) for Sx at given s. Due to the 
approximation in eq. Sj, s unless s = Seq- An 
approximation of Sx{s) is given in Table |lj The Y{s) ob¬ 
tained from the data of Fig. are shown in Fig. for 
different sizes. Symbols mark the results at s = Seq- 
Notice the strong increase of Y (s) for large sizes. For 
N = 37888 {Lq ~ 315 A), Y increases from ^ 100 N/m 
at Seq — 0.9985 to 220 N/m at s = 0.9995, i.e. more than 
a factor 2 for a strain e = s — Seq = 0.001 (0.1 %)! This 
strong dependence is in full qualitative agreement with 
the recent experimental claims [UlIIS]. A quantitative 
comparison will be discussed below. 

Since Beq and C'ii,eg, as well as Seq, C, D, C and D 
turn out to depend smoothly on Lq, we can approximate 
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FIG. 4. Young modulus Y(s) as a function of s for sample 
sizes as in Fig. 2 with symbols for s = Seq. The inset shows 
Y{s) for a size of 1 fj.m {N ~ 3.8 x 10^), calculated using eqs. 
|10| and and the expressions in Table [I] The upper axis of 
the inset gives the strain e = s — Ssq. 



FIG. 5. Young modulus Y as a function of Lq for the indi¬ 
cated values of the pressure P. The value in brackets is the 
corresponding strain e = s — Seq- 
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TABLE I. Size d epen den t p arameters for B{s) and Cii(sj;) 
according to eqs. |l0| and [l^ for Lq in A. Beq, Cn^eq, C and 
C are in eV/A^, other quantities are dimensionless. 


all parameters by the functions of Lq given in Table I. 
These expressions apply to any size and give appropriate 
asymptotics with C {C) equal to dP/ds {dPx/dsx) for 
the smallest system (iV=24). 

The inset of Fig. shows the resulting Y(s) for a size 
Lq = 1 /rm N = 3.8 x 10^ atoms). At zero strain 
(symbol) Y is only 30 N/m, becoming almost a factor 
10 larger at 0.5 % tensile strain, where it approaches 
its asymptotic value. Although the suppression of an- 
harmonicity found here goes very fast as a function of 
strain, it clearly deviates from membrane theory within 
the SCSA, where a complete suppression of anharmonic- 
ity occurs for tensile strain of 0.01 % [22]) two orders of 
magnitude lower than our atomistic simulations. 

Finally, the size dependence of Y with tensile strain at 
negative pressures is displayed in Fig. Tensile stress of 
0.05 N/m, corresponding to ~0.05 % strain, suppresses 
the anharmonic effects, and thus the power law decay, 
for Lq > 0.25 /im. As a consequence, Y is a factor ~ 4 
larger than Y^q for a system of 1 /im. Subsequently, in¬ 
creasing the stress by a factor 10 yields a strain of ~0.25 
% and Y ~ 275 N/m. This variation of Y with strain 
corresponds to recent experimental data. m The factor 2 
difference in both lower and upper bound of Y, however, 
with an experimental upper value Y =700 N/m, remains 


unexplained and requires further investigations. We note 
that our values of the bare elastic moduli agree with the 
experimental phonon spectrum 1301 ED of graphite, where 
anharmonic effects are suppressed by interlayer interac¬ 
tions. 

In conclusion, we have shown by atomistic simulations 
that the in-plane elastic moduli of graphene vanish with 
size as Lq^'^ with rju — 0.325, confirming that graphene 
follows membrane theory within the SCSA in this re¬ 
spect. The critical exponent ? 7 „, together with the inde¬ 
pendent estimate of ~ 0.85 jj, supports the scaling 
relation ?/„ = 2 — 2r]. In contrast, our results do not 
support the SCSA prediction i/ = —1/3 for L —b oo. We 
remark that v = —1/3 in eqs. [^andj^leads to Br = —Xr 
and Xr = 2Br — Ch, implying that Xr should be neg¬ 
ative for stability while in Fig. ^ 2Br — Cu^r remains 
positive for any Lq. We also find that suppression of an- 
harmonicity requires a tensile strain 10-50 times larger 
than predicted by SCSA. 
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